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We describe and analyze CUDA simulations of hydro dynamic interactions in active dumbbell 
suspensions. GPU-based parallel computing enables us not only to study the time-resolved collective 
dynamics of up to a several hundred active dumbbell swimmers but also to test the accuracy of 
effective time-averaged models. Our numerical results suggest that the stroke-averaged model yields 
a relatively accurate description down to distances of only a few times the dumbbell's length. This 
is remarkable in view of the fact that the stroke-averaged model is based on a far-field expansion. 
Thus, our analysis confirms that stroke-averaged far-field equations of motion may provide a useful 
starting point for the derivation of hydrodynamic field equations. 
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I. INTRODUCTION 

The derivation of effective hydrodynamic equations 
from microscopic or mesoscopic models presents a key 
problem of non-equilibrium statistical physics [T| . Stan- 
dard techniques typically involve severe approximations 
such as, for example, factorization of correlation func- 
tions, truncation of hierarchies, closure conditions, etc. 
Understanding the details of such approximations is cru- 
cial for identifying the range of applicability of the re- 
sulting field equations. If a complex fluid is made up 
of active constituents (e.g., bacteria or other micro- 
organisms) that propel themselves by quasi-periodic 
swimming mechanisms [2J [3J, then one usually faces 
the additional task of approximating the explicitly time- 
dependent microscopic dynamics with a set of coarse- 
grained, time-averaged equations of motion. Aiming at a 
better quantitative understanding of this commonly em- 
ployed approximation, the present paper provides a de- 
tailed comparison between the microscopic dynamics of 
actively driven, spring-based dumbbells and those of a 
time-averaged analytic model derived from far-field ex- 
pansion [H[S]. 

Owing to the fact that hydrodynamic interactions are 
long-range, simulations of the full time-resolved dynam- 
ics of S = N/2 dumbbells (each consisting of two spheres) 
are numerically expensive, scaling as TV 2 . However, in the 
deterministic limit case and/or additive noise limit [5J, 
the dynamics is well suited to parallel computations. 
Very recently, GPU-based codes have been used for var- 
ious statistical mechanics problems, yielding speed-ups 
on the order of 20-100 times a CPU-only solution [TUTU]. 
Here, we implement a straightforward N 2 solution of the 
hydrodynamic equations of motion, a communications- 
intensive task which is difficult to parallelize in tradi- 
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tional clusters. For a moderate population size (up to a 
few thousand particles), this method decreases the com- 
putation time by a factor of 100 compared with conven- 
tional CPU simulations on standard consumer hardware. 
Hence, we identify GPU computing as a promising ap- 
proach for future simulations of active particle suspen- 
sions (details of the numerical implementation are sum- 
marized in Sec. jv|. 

Passive (non-driven) dumbbell models have been 
widely investigated in polymer science and related fields 
over the past decades (see, e.g., Refs. [ITtiTT] ). Very re- 
cently, several authors [U [TS] considered active, inter- 
nally driven dumbbells as prototype systems for collec- 
tive swimming at zero Reynolds number, 51 = 0. Loosely 
speaking, one can say that active dumbbell systems con- 
stitute a sort of 'Ising model' of collective swimming, i.e., 
they represent strongly simplified models which can be 
treated by analytical means, thus providing useful in- 
sights. Active dumbbells are particularly well-suited to 
identifying the role of hydrodynamic long-range interac- 
tions in collective micro-swimming. This is because iso- 
lated deterministic dumbbells are prevented from self- 
motility by Purcell's scallop theorem [3J. Hence, any 
effective motion in deterministic dumbbell systems is 
caused by hydrodynamic interactions between different 
dumbbells. 

In a recent paper, Alexander and Yeomans [5] have de- 
rived analytical expressions for the effective far-field in- 
teractions of symmetric, active dumbbells in three dimen- 
sions (3d). Specifically, they showed that the effective 
hydrodynamic pair interaction decays with distance \D\ 
as |-D| -4 - Considering Id motions, Lauga and Bartolo [4 
extended this result to asymmetric dumbbells and found 
that in this case the hydrodynamic interaction decays 
less strongly as |-D| -3 - While these studies have led to 
novel insights into interplay between swimmer symmetry 
and effective long-distance interaction scaling, a detailed 
comparison of microscopic and stroke-averaged models is 
still lacking. The present paper intends to close this gap 
with respect to symmetric dumbbells. 

For this purpose, we shall first introduce a microscopic 
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spring-based dumbbell model (Sec. [IT]) that can be readily 
implemented in GPU-based computer simulations. In the 
limit of an infinitely stiff spring, our model reduces to a 
shape-driven dumbbell model as considered in Refs. |H[S]. 
The corresponding 3d stroke-averaged equations of mo- 
tion will be discussed in Sec. |III| After having confirmed 
that the stroke-averaged model reproduces the main fea- 
tures of the microscopic model simulations in Id, we per- 
form similar comparative studies for 3d arrays of sym- 
metric dumbbells. Generally, we find that the stroke- 
averaged dynamics yields relatively accurate description 
of the microscopic model down to distances of only a 
few times the dumbbells' length. This is remarkable in 
view of the fact that the stroke-averaged model is based 
on a far-field expansion. Thus, at least for the model 
considered here, our results suggest that stroke-averaged 
far-field interaction models may indeed provide a useful 
starting point for the derivation of hydrodynamic field 
equations. 



II. MICROSCOPIC MODELING OF ACTIVE 
DUMBBELLS 

We shall begin by summarizing the "microscopic" 
model equations of the spring-based dumbbells simulated 
in our computer experiments. The corresponding stroke- 
averaged equations of motion will be discussed in Sec. |III| 
To keep the discussion in this part as general as possible 
- and as reference for future work - we shall formulate 
the model for "Brownian" dumbbells, even though the 
discussion in the subsequent sections will be restricted to 
the deterministic limit. 

We consider a system of S identical dumbbells. Each 
dumbbell consists of two spheres, of radius a. At very 
low Reynolds numbers, inertia is negligible and the state 
of the system at time t is completely described by the 
spheres' position coordinates {-X'a} = {X/ ai -\(t)} with 
a = 1, . . . , 2S labeling the spheres, and i = 1, 2, 3 the 
space dimension (throughout, we adopt the Einstein sum- 
mation convention for repeated lower Latin indices) . Ne- 
glecting rotations of the spheres, the dynamics is gov- 
erned by the Ito-Langevin equations [HI HM23] 
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where ks denotes the Boltzmann constant and T the tem- 
perature of the surrounding fluid (X :— dX/dt). Equa- 
tion (la I corresponds to the overdamped limit of Stoke- 
sian dynamics [23]. The Gaussian white noise £( 7 fc)(t) 
models stochastic interactions with the surrounding liq- 
uid molecules and is characterized by |25j 
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The hydrodynamic interaction tensor !K couples the de- 
terministic force components -F^i) that act on the indi- 
vidual spheres. Generally, the vector F — {F^} com- 
prises contributions from internal forces, e.g., those re- 
quired to bind and oscillate spheres in an active dumb- 
bell, as well as from external force fields (gravity, etc.). 

The amplitude of the noise force is determined by the 
fluctuation dissipation theorem, which is satisfied if C is 
constructed from "K by Cholesky decomposition, i.e., 
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In our numerical simulations, "K is given by the Rotne- 
Prager-Yamakawa-Mazur tensor 
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where r a pi := x ai - xpu a ^ f3, and r af} := \x a - xp\. 
Analytical formulas presented below are based on an 
Oseen approximation, which neglects the second line in 
Eq. (3b). The diagonal components (3a) describe Stoke- 



sian friction in a fluid of viscosity /x. The off-diagonal 
components ( 3b I model hydrodynamic interactions be- 
tween different spheres. Note that "K is positive defi- 
nite for r Q( g > 2a and divergence-free, d(pj)yt(ai)(Pj) = 
with <9(/3i) := dldxtpi), implying that the Cholesky- 
decomposition ([2| is well-defined. 

To completely specify the model, we still need to fix 
the intra-dumbbell force. To this end, consider the dumb- 
bell a, formed by spheres a = 2a — 1 and (3 — 2a, and 
denote its length by d a (t) := \Xp(t) - X a (t)\. Neglect- 
ing external force fields from now on, we shall assume 
that the two spheres are connected by a harmonic spring 



of variable length L a (t), i.e., F( 
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d(/3i)U where 
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with L a {t) — I + \ sm(uit + ip 17 ) denoting the time- 
dependent equilibrium length of the spring, and £ the 
mean length such that £ > 2a + A. The dumbbell swim- 
mer is called passive if the stroke amplitude A = 0, and 
active if |A| > 0. As discussed below, the phase param- 
eter ip a is important for the interaction between two or 
more dumbbells. 

For the overdamped description ([!]) to remain valid, 
the driving must be sufficiently slow. More precisely, we 
have to impose that T 7 <C To <C T, where T := 2tt/lu 
is the driving period, To := 2n j \Jk^jM the oscillator 
period for a sphere of mass M, and T 1 := M/7 the char- 
acteristic damping time. This restriction ensures that 
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the dumbbells behave similar to shape-driven swimmers, 
i.e., d? ~ L° (t) is a useful approximation in analytical 
calculations. 

With the above assumptions, the TV-particle PDF 
f(t,{xr a i)}) of the stochastic process {-X7 ai )(t)} from 
Eq. ([!]) is governed by the Kirkwood-Smoluchowski equa- 
tion 
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Here T = 2it/uj denotes the period of a swimming stroke. 
By assuming that TV (t) and R (t) are slowly varying 
functions of time, one can further approximate 
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for any sufficiently well-behaved function /. 



For time-independent potentials, the stationary solution 
of this equation is given by the Boltzmann distribu- 
tion, / oc e -u/(k B T) _ However, in the remainder, we 
shall focus on the deterministic limit case, formally ob- 
tained by putting T = in Eqs. (la I, which is justified 
for sufficiently big spheres. 



III. STROKE-AVERAGED HYDRODYNAMIC 
PAIR INTERACTIONS 

In this part we will summarize the stroke-averaged 
equations of motion for the case of 3d symmetric, de- 
terministic dumbbell swimmers (a detailed derivation, 
which differs slightly from that in Ref. [3] but yields 



equivalent results, is given in the Appendix). In Sec. IV 



the dynamics resulting from these effective equations of 
motion for the dumbbell positions and orientations will 
be compared with numerical simulations of the micro- 
scopic model equations ([I]). From now on all considera- 
tion refers to the deterministic limit case. 



A. General stroke-averaging procedure 

We characterize each dumbbell by its direction vector 

X 2a — X 2a -1 



N (t) 
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and its geometric center 
1 



(7=1, 



R (t) := -(X 2a + X 2a _ l ). 
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Note that for symmetric dumbbells the geometric center 
coincides with the center of hydrodynamic stress [TSJ [2j5] . 

The basic idea of the stroke-averaging procedure [U 
5, 23 is to focus on the dynamics of averaged position 
and orientation coordinates R(t) and N a (t), which are 
defined by 

N°(t) := - / duN (u), (7a) 

1 Jt-T/2 
nt+T/2 

R°(t) := - / Auil(u). (7b) 

1 Jt-T/2 



B. Stroke-averaged equations of motion 

Using the approximations (18b, one can derive from the 
microscopic model equations (Tib with T = the corre- 
sponding deterministic, stroke-averaged, far-field equa- 
tions of motion [3] [3J [53] , by making the following sim- 
plifying assumptions: 

(i) The dumbbells are force-free and torque- 
free|41j and approximately shape-driven, i.e., 
d a := \X 2 a~X 2a _ 1 \ ~e + \sm(ut + <p a ). 

(ii) The dumbbells are slender, i.e., sphere radius a 
and stroke amplitude A have about the same size, 
but are much smaller than the dumbbell's mean 
length I. 

(iii) The ensemble is dilute, meaning that the distance 
D a P ._ j^Pj -=\r° - rp\ between dumbbells a 
and p is much larger than £. 

Adopting (i)-(iii) and restricting to two-body interac- 
tions, one finds the effective equations of motion 
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where the stroke-averaged hydrodynamic interaction 
terms to leading order in \jt are given by 
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Here, the unit vector D P := D ap /\D ap \ gives the orien- 
tation of the distance vector D ap = R a — R p , and s, r, q 
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abbreviate the projections 



— n a P 



D] P N?, q := NJN? 



(9d) 



One readily observes the following prominent fea- 
tures: (i) The effective translational interactions scale 
as oc \D ap \~ i . (ii) The effective rotational interactions 
scale as oc \D ap \~ 5 . (iii) The stroke-averaged interaction 
terms J, K vanish if the phases <p a and ip p differ by mul- 
tiples of 7r [5]. This illustrates the importance of phase 
(de)tuning in the collective swimming at zero Reynolds 
number. 



IV. MICROSCOPIC VS. STROKE-AVERAGED 
DYNAMICS 

We next compare the predictions of the stroke- 
averaged equations ([9| with numerical results obtained 
from CUDA simulations of the microscopic spring-based 
dumbbell model from Sec. [IT] For this purpose, we first 
consider Id aligned dumbbell pairs similar to those stud- 
ied by Lauga and Bartolo [4]. The rotational interaction 
of two dumbbells will be analyzed in Sec. |IVB| Finally, 
we also study the collective motion of 3d grids of dumb- 
bells (Sec. IV C). In all cases, the swimmers are assumed 
to be in an infinite body of fluid initially at rest, i.e., no 
additional boundary conditions (periodic or otherwise) 
are imposed. 



The quantity R 21 (t) characterizes the net motion of the 
dumbbell pair, whereas AR 21 (t) indicates whether the 
dumbbells the move towards or away from each other. 

As is evident from Fig. [I] (a), symmetric dumbbells 
move in the same direction with identical speeds; the 
direction of the motion is determined by the phase dif- 
ference ip 2 — ip 1 . As predicted by Eq. (10), the col- 



lective displacement over a swimming stroke varies as 
|Z? <Tp |~ 4 with the distance between the dumbbells, see 
Figure |T| (b) . Even though the stroke-averaged equa- 
tions ( 10 1 are based on a far-field expansion, they de- 



scribe the microscopic dynamics of aligned dumbbells 
well down to distances of a few body lengths. 

In this context, it is worthwhile to note that the quan- 
titative difference between the stroke-averaged dynamics 
(solid lines) and the microscopic simulations (symbols) 
in Figs. [I] (a) and (b), is due to the relatively large pa- 
rameter ratio a/£ ~ 0.2 used in these simulations. As 
shown explicitly in the Appendix, the stroke averaged 
equations of motion ( |10| become more accurate in the 
limit a/£ — > 0. This is confirmed by the numerical re- 
sults shown in Fig. [T (c) . This diagram depicts the ra- 
tio of the average co lective swimming speeds (i.e., the 
collective displacements after a stroke period) obtained 
by either method for different values of a/i at constant 
stroke-amplitude A. We readily observe that this ratio 
approaches unity in the limit a/£ — ¥ 0. However, in view 
of the fact that the collective swimming speed is approx- 
imately proportional to the sphere radius a, see Eq. (10), 



A. Aligned dumbbell pairs 



we opted for a moderate value ajl = 0.2 in all our simu- 
lations in order to observe noticeable swimming effects. 



As long as thermal fluctuations are negligible, aligned 



dumbbells do not change their orientation and Eq. (9a 
reduces to (see A 1 for an explicit derivation) 
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where R a denotes the coordinates along the common 
axis. The lines in Figs, [I] (a) and (b) represent the dy- 
namics of aligned dumbbell pairs (S = 2) as predicted 
by Eq. (10). Symbols were obtained from microscopic 



simulations of the corresponding spring-based model de- 
scribed in Sec. [TlJ Following Lauga and Bartolo [I], we 
quantify collective motion of the dumbbell pairs in terms 
of their mean collective displacement (solid lines/filled 
symbols in Fig. [T]) , 



RK(t) = -[R 2 (t) 



■R\t)], 



(11a) 



and their mean relative distance (dashed lines/unfilled 
symbols in Fig. [T]) , 



AR 21 (t) = R 2 (t)-R 1 (t). 



(lib) 



B. Two-dimensional rotation of dumbbell pairs 

Dumbbells that are arranged in an aligned Id config- 
uration do not change their orientation. This is differ- 
ent for non-aligned configurations in higher dimensions 
where hydrodynamic pair interactions can induce rota- 
tions. To test the accuracy of the stroke-averaged equa- 
tion (9b) for the rotational motions in two dimensions, 



we conducted a series of simulations using the following 
setup: The first dumbbell (labelled by a) was placed at 
the origin oriented along the x-axis, and a second dumb- 
bell (p) was placed such that the geometric centres were 
separated by a distance of hi. By varying the starting po- 
sition of the second dumbbell along a circle, while keeping 
the initial projection constant, we can compare numeric 
and analytic results for various projections s(t), r(i), q(t), 
as defined in Eq. (9d). 



Figure [2] depicts the change of the dumbbells' relative 
orientation 

Aq(t) := q(t) - g(0) , q(t) := N?(t)N?(t) (12) 

after five swimming strokes t = 5T for two different ini- 
tial projections (a) q(0) = and (b) q(Q) — 1. As evident 
from the diagrams, in both cases the stroke averaged de- 
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FIG. 1: Comparison of exact microscopic motion and effective stroke-averaged dynamics for an aligned dumbbell pair. Lines 
were obtained by numerical integration of the stroke-averaged equation |l0f , whereas symbols show the simulation results for 
the microscopic spring-based dumbbell model described in Sec. [IT] Solid lines and filled symbols depict the mean displacement 
B^it) - fl2T(0) = {[R 2 (t) + R 1 ^)} - [R 2 (0) + R 1 (0)]}/2 of the geometric centres. Dashed lines and unfilled symbols indicate 
the relative separation AR 21 (t) - AR 21 (0) = [R 2 {t) - - [i? 2 (0) - ^(O)]. (a) Symmetric dumbbells do not change their 

relative separation and move linearly in time depending on the phase difference Aip = ip 2 — tp 1 . Simulation parameters are 
comparable to those of Lauga and Bartolo [3J: Initial separation Ai? 21 (0) = -R 2 (0) — ^? 1 (0) = W£, mean dumbbell length 
t — 5/im, driving frequency lj — 500s _1 (time on the x-axis is given in units of the stroke period T = 2-k/u), stroke amplitude 
A = 0.1£, a = 0.2£, phase difference tp 2 -— tp 1 = tt/2. For the microscopic model: spring constants fco = 0.001kg/s 2 , viscosity 
/i = 10~ 3 kg/(ms), particle mass density g — 10 3 kg/m 3 , simulation time step At « 10 _4 T. (b) Distance dependance of the 
collective motion and separation for aligned dumbbell pairs during a stroke period T. Line styles and symbols corresp ond 
to the same configurations and simulation parameters as used in (a). Remarkably, the stroke-averaged far- field equation (10 1 
describes the microscopic dumbbell dynamics well down to distances of a few body lengths; however, the deviations from the 
time-resolved microscopic dynamics accumulate over time, as is evident from (a). The difference between the stroke-averaged 
dynamics (solid lines) and the microscopic simulations (symbols) in (a) and (b) is due to the choice of a relatively large 
parameter ratio a/l = 0.2 in these simulations; the results of both methods agree in the limit a/l — > as illustrated in diagram 
(c), which shows the ratio of mean swimmer displacements obtained from the microscopic ('Mic') and stroke-averaged ('SA') 
dynamics at constant A = 0.1^ and various choices of a/l. 



scription ( 9b | correctly reproduces the rotational dynam- 
ics of the microscopic spring-based model. 

We may thus briefly summarize: The results in Figs. [T] 
and [2] show that the stroke-averaged equations ([£]) sat- 
isfactorily capture the main features of effective pair 
interactions in the spring-based microscopic model at 
moderate-to-low densities (large distances). This cor- 
roborates that equations of the type (|9) can provide a 
convenient mesoscopic description which, for example, 
can be used as a starting point for derivation of coarse- 
grained macroscopic field theories [31]. Conversely, the 
good agreement between the averaged dynamics ([9]) and 
the microscopic model simulation provides a helpful con- 
firmation that our CUDA algorithm works correctly even 
at relatively low densities, when hydrodynamic interac- 
tions effects are relatively weak and algorithms may be- 
come prone to numerical instabilities. 

In the remainder, we shall focus on 3d many-swimmer 
simulations that fully exploit the virtues of the CUDA 
parallelization scheme. 



C. Collective swimming of three-dimensional 
dumbbell arrays 

In this section we will compare the predictions of the 
stroke-averaged far-field equations ^ with simulations of 
spring-based dumbbells for 3d dumbbell configurations. 



In our simulation the dumbbells' geometric centers 
R a (0) are initially placed on a cubic (x x x x x)-latticc 
with equidistant spacing p -1 / 3 , where p is the number 
density of the configuration. Initial orientations N a (0) 
are sampled uniformly from the unit sphere. For the lat- 
tice size we consider values x = 3, 5, 7, 9 corresponding to 
a total dumbbell number S = 3 3 ,5 3 ,7 3 ,9 3 , respectively. 
To characterize the collective motion, we measure in our 
simulations the mean square displacement per particle 
averaged over different initial conditions, i.e., 



w s 

^) 2 » := w E s £ [Ra{t] w) ~ Ra{0; ' 
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(13) 



where the variable w = 1, . . . , W labels different initial 
conditions. We distinguish two classes of initial condi- 
tions: 

(i) An "optimized" phase distribution: Phases were 
set such that each dumbbell had a phase of or 
7r/2 with all nearest neighbors having the alternate 
phase in the manner of a 3d "checkerboard" . The 
corresponding results for the microscopic simula- 
tion and the stroke-averaged model are indicated 
by filled symbols and solid lines in Figs. [3] and |4j 
respectively. 
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FIG. 2: Rotational motion of dumbbell pairs; symbols indi- 
cate numeric data while lines represent analytics. 6^ is the 
angle between a line connecting the dumbbell's geometric cen- 
tres and the :r-axis; Aq(5T) := q(5T) — q(0) where q(t) is the 
projection of the swimmer orientations q(t) = Nj(t)Nj(t). 
Initial configurations: (a) q(0) = and (b) q(0) = 1 with 
an initial radial separation of 51 . One readily observes the 
good agreement between the microscopic simulations and the 
analytics. 

(ii) A randomized phase distribution: Phases were set 
to random values evenly distributed on the interval 
[0, 2tt) , with a different distribution for each run. 
The corresponding results are indicated by unfilled 
symbols and dashed lines in Figs. [3] and |4j respec- 
tively. 

Fig. [3] illustrates how the mean square displacement 
over a period, ((R(T) 2 )}, varies with density p - or equiv- 
alcntly with grid spacing - for an array of (5 x 5 x 5) = 125 
dumbbells. As evident from the diagram, the predic- 
tion from the stroke-averaged model (|9| is in good agree- 
ment with the scaling behavior measured for the micro- 
scopic model. Furthermore, by comparing filled with un- 
filled symbols and solid with dashed lines, we note that 
the collective displacement ((R(T) 2 )) is generally smaller 
for the randomized phase distribution, corroborating the 
fact that optimizing the phase distributions can consider- 
ably enhance the effectiveness of collective motions [35] • 

Figure [I] shows how the quantity ((R(T) 2 )) scales with 
the total number S of the dumbbells at fixed density p. 
After a slight initial jump from the (3x3x3) case, 
adding more swimmers at fixed density p produces only a 




FIG. 3: (a) Scaling of the mean square displacement per par- 
ticle (measured over a period) with number density for two 
different phase distributions. The diagram depicts the simu- 
lation results for a cubic array of (5 x 5 x 5) = 125 dumbbells, 
averaged over W = 100 different runs, each with random ini- 
tial orientations. Symbols refer to the spring-based model and 
lines to the stroke-averaged model J9J. The collective mean 
square displacement is proportional to {ip 1 ^ 3 ) 2 " with an ex- 
ponent v — —4. (b) Mean square displacement rescaled (i.e., 
multiplied) by p~ 8//3 . We observe that for an "optimized" 
phase distribution (filled symbols/solid lines) the effective 
mean square displacement is larger than for a uniformly ran- 
dom phase distribution (empty symbols/dashed lines). On 
this scale, statistical error bars (not shown) are smaller than 
the size of the symbols. The shift between lines and symbols, 
caused by the relatively large parameter ratio a/ 1 = 0.2 used 
in these simulations, is consistent with the value expected 
from Fig. [T] (c). 



minimal increase in displacement, and the effect appears 
for both optimized or randomized phase distributions. 
Again, collective displacement is generally smaller for 
randomized phase distributions than for optimal phase 
distributions. 
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FIG. 4: Scaling of the mean square displacement per parti- 
cle over a period with dumbbell number S at constant den- 
sity. The diagram depicts the simulation results for collec- 
tions of dumbbells arranged on a cubic (x x x x a^)-lattice 
with x = 3, 5, 7, 9 and spacing 101, averaged over W = 10 
different random initial orientations. Symbols refer to the 
spring-based model and lines to the far-field stroke-averaged 
model (pll, beginning from the same initial conditions. In- 
creasing the number of swimmers while keeping the number 
density constant produces only minimal gains in translational 
speed. The collective mean square displacement is smallest 
for (3x3x3) = 27 swimmers, which is due to the rela- 
tively large number of swimmers with an incomplete set of 
"nearest neighbors" in this case. We also observe that for 
an "optimized" phase distribution (filled symbols/solid lines) 
the effective mean square displacement is larger than for a 
uniformly random phase distribution (empty symbols/dashed 
lines). Again, the shift between lines and symbols, caused 
by the relatively large parameter ratio a/l = 0.2 used in 
these simulations, is consistent with the value expected from 
Fig- 0(c). 



V. COMPUTATIONAL ASPECTS 

The numerical results were obtained from parallelized 
simulations run on graphics processing units (GPUs) 
using Nvidia's Compute Unified Device Architecture 
(CUDA). Compared to conventional CPU programs, 
GPU algorithms may yield significant speed ups (up 
to factors of a few hundreds) whenever a problem can 
be naturally parallelized [TMlO]. on relatively low-cost 
consumer-grade hardware. In cases where the problem is 
small enough to fit on a single device, the resulting soft- 
ware is simpler, easier to test, less costly to implement, 
and much faster than traditional cluster-based methods. 
This is the case for deterministic many-swimmer simula- 
tions, for stochastic single-swimmer simulations, and also 
for stochastic many-swimmer simulations with purely ad- 
ditive noise, corresponding to a constant matrix C in 
Eq. 0. 

Most 0(N 2 ) problems such as TV-body simulations 
with pair interactions involve enough data communica- 
tion that they are difficult to distribute efficiently, or 



are costly enough that they must be recast in more nu- 
merically tractable forms such as Ewald summation (33j- 
155] . For cases of a few hundred swimmers, CUDA- 
based implementations of straightforward 0(JV 2 ) prob- 
lems present an excellent method for numerical simu- 
lation. We use a simple direct computation of sphere- 
sphere interactions via the Rotne-Prager-Yamakawa- 
Mazur tensor, disregarding lubrication forces and close- 
range interactions due to the dilute nature of the suspen- 
sions and slender structure of the dumbbells. 

We tested our GPU code on an AMD Phenom X4 
940 system running Fedora Linux, using a consumer-level 
Nvidia GTX 295 GPU and a more research-oriented Tesla 
C1060 GPU; other tests took place on an Intel i7 860 run- 
ning Gentoo Linux and a consumer- level GTX 276 GPU. 
All hardware was capable of double-precision calculation 
and used version 2.3 of the CUDA toolkit. 

Initial testing of a similar but simpler problem (colloids 
moving under constant applied force, using Oseen inter- 
actions and single precision) showed very large bench- 
marked speed-ups compared with a C-based CPU simu- 
lation. For example, with ~ 2000 particles, we measured 
up to a ~ 450 x speed-up when calculating the full hy- 
drodynamic interaction tensor and ~ 800 x speed-up us- 
ing an un-optimizcd version of the elegant tiled method 
described in We did not benchmark the current 

simulation, but estimate the speed-up, while still being 
significant, to be considerably less due to the complexity 
of the multi-swimmer problem, additional data transfers 
from the GPU, and the use of double precision. 

Despite the speed advantage of the tiled method, we 
decided to compute the full hydrodynamic interaction 
tensor in our simulations, primarily to maintain congru- 
ence with existing C code in a battery of automated unit 
tests and to keep the code as simple as possible. Future 
implementations will likely reintroduce the tiled calcula- 
tion to further improve computational efficiency. Gener- 
ally, it is encouraging that even a relatively straightfor- 
ward CUDA simulation of the multi-swimmer problem 
exhibits compelling speed advantages over a CPU-based 
solution. 

The use of double precision is unfortunate in that single 
precision calculations on CUDA processors show signifi- 
cant performance increases due to better hardware sup- 
port and memory performance. However, in the case of 
collective dumbbell motion, the distance moved in each 
time-step is very small compared to the length of the 
dumbbells or their position, which caused initial calcula- 
tions using single precision to fail, as the position incre- 
mental during a single timestep fell below the threshold of 
machine precision. To allow for standardized testing, we 
elected to use double precision and to accept decreased 
performance rather than implementing a better accumu- 
lation algorithm (such as Kahan summation) based on 
single precision. For reasons of accuracy, we also chose 
not to enable Nvidia's fast-math optimizations. The lat- 
ter can significantly accelerate the computation of certain 
numeric functions (particularly trigonometric functions) 



but this gain comes at the cost of some precision. How- 
ever, this might represent another opportunity for per- 
formance optimization in the future. 

Another important issue is the choice of the integra- 
tor due to accumulation of errors and the stiffness of the 
problem. A variety of methods were tested, including 
Euler, Adams-Bashforth-Moulton, and Runge-Kutta in- 
tegrators. The approach eventually used was a one-step 
Heun predictor-corrector method, which produced ex- 
cellent results and can easily incorporate additive noise 
for stochastic simulations. The time step for the sim- 
ulations was chosen based on the smallest dynamical 
time scale in the problem (given by the spring frequency 
To = 2n/^/M, see discussion in Sec. Ilh and then 
manually reduced until numerical errors were acceptable 
by ensuring that single dumbbells did not translate and 
numeric fluctuations were orders of magnitude below the 
expected motion caused due to hydrodynamic interac- 
tions. 

The use of a spring-based model created an addi- 
tional complication: After prescribing the initial posi- 
tion, orientation, and phase of the dumbbell we initially 
placed the spheres centred at the potential minima. How- 
ever, numerical integration and finite potential strengths 
caused the sphere positions to lag very slightly behind the 
potentials once they began moving periodically. Since 
the hydrodynamically induced dumbbell motion is of a 
very small scale compared to the dumbbell size, this ini- 
tial settling caused a large anomalous motion during the 
first period of simulation. To rectify this, it was neces- 
sary to discard the first period and begin measurements 
after the lag was established and dumbbell translation 
was approximately linear. This did cause miniscule de- 
viations of the dumbbells' mean length £, amplitude A, 
and phase ip from the values specified by the initial con- 
ditions, but tuning the potential spring constant to be 
sufficiently stiff reduced these deviations to acceptable 
values of a few percent. 

While the results shown here are purely deterministic, 
incorporating noise is relatively straightforward, as the 
system hydrodynamic tensor 5£ may be numerically de- 
composed via Cholesky decomposition [37l [38] . However, 
even with GPU acceleration this decomposition is pro- 
hibitively expensive; in the case of slender dumbbells and 
dilute suspensions, we advocate a simple additive noise 
with a constant matrix C as a reasonable approximation 
in the dilute limit as the off-diagonal terms in Eq. ^ are 
negligible. We compared the full Cholesky decomposi- 
tion and an additive-noise approximation in various test 
runs and found that the results for the collective mean 
square displacement differed by only a few percent. 



VI. CONCLUSIONS 

We have examined the stroke-averaged, far-field equa- 
tions of motion for symmetric dumbbells, and veri- 
fied the general properties of this coarse-grained model 



by comparing with microscopic numerical simulations 
at relatively low densities. Remarkably, the micro- 
scopic and coarse-grained simulations agree well even at 
intermediate-to-high swimmer densities, where the effec- 
tive equations of motion are expected to become less ac- 
curate. However, it should be kept in mind that at very 
high densities, when collisions (i.e., steric effects) become 
relevant, lubrication effects as well as near-field hydrody- 
namics must be modelled more carefully. 

In the case of dumbbells arranged on a 3d grid, the 
translational speed due to hydrodynamic interaction be- 
tween dumbbells varies predictably with spacing, tend- 
ing toward |-D| -4 decay, where \D\ is the distance be- 
tween dumbbell centres. Due the short range of the ef- 
fective hydrodynamic interactions for symmetric dumb- 
bells, adding more swimmers at a fixed density has only a 
minimal impact on dumbbell translational speed. On the 
other hand, the collective swimming speed can be notice- 
ably increased by replacing a randomized phase distribu- 
tion with an ordered, "optimal" distribution of phases 
such that the difference in phase between a periodically- 
driven dumbbell and its nearest neighbors is tt/2. 

Generally, our numerical investigations illustrate that 
GPU-bascd simulations of multi-swimmer systems can 
provide a valuable tool for studying collective motions at 
very low Reynolds number. Moreover, the CUDA algo- 
rithm used in our computer experiments can be readily 
adapted to simulate hydrodynamic interactions between 
colloids that can be trapped and manipulated by means 
of optical tweezers [35] • Such theoretical investigations 
can help to create more efficient micropumps, e.g., by op- 
timizing the phase relations in oscillating arrays of col- 
loids. 

Finally, another long-term objective is to compare 
many-swimmer simulations with predictions of effective 
field theories [10]. Our above results suggest that the 
most promising approach towards achieving this goal 
may be a two-step procedure: (Step 1) One should try to 
derive stroke-averaged equations of motion that correctly 
capture the phase dependence on the level of effective 
two-particle interactions. As our above discussion has 
shown, such coarse-grained models can correctly repro- 
duce many of the main features of the microscopic model. 
Thus, it is sufficient for many purposes to implement the 
coarse-grained equations into a CUDA environment (step 
2). Compared to simulations of the full microscopic dy- 
namics, this may reduce the effective simulation time by 
an additional factor of 100 or more since the analytic 
stroke-averaging procedure makes it unnecessary to nu- 
merically resolve the smallest dynamical time scales in 
the system. We hope that our analysis may provide use- 
ful guidance for future efforts in this direction. 
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Appendix A: Stroke-averaging 

We derive the stroke-averaged effective interactions 
between two symmetric, quasi-shape-driven dumbbell 
swimmers. Each dumbbell consists of two spheres of ra- 
dius a. The swimming stroke of an individual dumbbell 
is assumed to be both force-free and torque-free. 



1. One-dimensional case 

In one space dimension (Id) we denote the position of 
the spheres belonging to dumbbell a by Xf , s = 1,2. To 
characterize position and orientation of the dumbbell, we 
may introduce center-of-mass and relative coordinates by 



K> = \{X<[+X%) (Ala) 

S a = X% - X{ (Alb) 
A CT = {Xl-Xl)l\X° 2 -Xl\- (Ale) 



hence, 



X{ = R" - S a /2 



XI = R a + S a /2 (Aid) 



which may also be written as 

X a s = R a + (-l) s S a /2. 



(Ale) 



Furthermore, we define the vector connecting two swim- 
mers a and p by 

D ap := R° -Rf (A2a) 
D ap := (R a - R p )/\R a - R p \. (A2b) 

The force-free constraint for the dumbbell a can be writ- 
ten as 



n = -n =■■ r, 



(A3) 



with F" denoting the internal forces acting on the first 
and the second sphere during a swimming stroke. Ne- 
glecting thermal fluctuations, the Id equations of motion 
can be written as 



(A4) 



p.r 



Here, we sum over all swimmers p = 1, . . . , S and the 
spheres r = 1,2 of each swimmer. Our goal is to de- 
rive from Eq. (A4) a stroke-averaged effective equation 
of motion for R" (for shape-driven dumbbells the motion 
of S a is trivial in Id). 



The "diagonal" components of the hydrodynamic in- 
teraction tensor H are given by the inverse Stokes friction 
coefficient 



i/ s 7= 7 - 1 = (67T Ala )- 1 . 



(A5) 



Adopting the Oseen approximation, the "off-diagonal'' 
components (s ^ r) read 



\S a \ 1 sr \X° ~x p \ 
where k = {Airp)^ 1 . It is useful to rewrite 

x°~ - x p r = D ap + r s 7 

where 



(A6) 



(A7a) 



(A7b) 



Using the force free condition (A3), we obtain from 
Eq. ([All 



(A8a) 



and 



S° 



p 

J2 B ap f p (A8b) 



Considering approximately shape-driven dumbbells, we 
have 



S a = L a (t)N a , 
\S°\ = L«(t), 
S a = L a (t)N a , 



(A9) 



where the periodic function L a {t) > describes the 
shape (length) of the dumbbell at time t. Hence, in- 
verting (A8b| we obtain the force as a function of the 
shape 



(A10) 



where B 1 denotes the inverse of the (S x S')-matrix 
B := (B ap ) defined in (|A8b|). Substituting this result 



into Eq. (A8a| yields the following closed equations for 
the position coordinates 



R a = J2 A ap (B- 1 ) pu L l/ (t) N v . 



(All) 



By means of Eqs. (A9), we can rewrite the off-diagonal 
components of the Oseen tensor as 

H™ = ^r, HZ P = K (Al2a) 



\D°p + Y° 
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where 



Y s 7 = -[(-l) s L°(t) N° - (-iy L p (t) NP] 

(Al2b) 

For a system consisting of more than two dumbbells 
(S > 2), the rhs. of Eq. (A10) contains not only two- 
body, but also three-body, four-body, . . . , S-body con- 
tributions. However, focussing only on the dominant 
two-body contributions, B := (B ap ) can be exactly in- 
verted and the rhs. of Eq. (All I can be expanded in 
the low-density limit corresponding to \D ap \ — > oo. Av- 
eraging the resulting power series over a stroke period 
[t - T/2, t + T/2] as described in Sec. |III A| and keeping 
only the leading order contribution, we find the follow- 
ing Id stroke-averaged equation of motion in two-body 
approximation 
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— ( 
16 



D°p\ 



D° 



2. Three-dimensional case 

In the 3d case, the derivation of stroke-averaged equa- 
tions becomes more complicated due to the additional 
rotational degrees of freedom. 

As before, we consider a dilute suspension of a = 
1,. . . ,N geometrically identical dumbbells of prescribed 
length L"{t). To characterize the motion of the dumb- 
bells, we define position and orientation vectors by 



R°{t) := \{X° l 



X a2 ) 



(A13) 



with S a denoting the non-normalized orientation vector, 
i.e., for a shape-driven dumbbell we can write 



S a (t) := X a2 - X c 
S a (t) := \S a \ = L a 



L"(t)N a , 



(A14a) 
(A14b) 



Similar to Eq. ( Ale ) , we can recover the bead coordinates 
{X al ,X a2 } from"{i2' T ) iV <T } by means of 

x *s _ R a + (_iy N ° L °/2. (A14c) 
As before, we consider shape-driven dumbbells with 



L a (t) = I + Asin(wi + tp a ). From the definition ( |A13[ >, 
one then finds that the exact equations of motion for 
{H CT ,iV' T } are given by 



1 rr(CTs)(pr) jppr 

s,p,r 



(A15a) 



N? = (S ik - N°N%) x 



(a2) (pr) 
kj 



-H 



(*l)(pr) 
kj 



p^a,r 



F pr 
J 



(Al5b) 



The indices s,r € {1,2} label the spheres and, through- 
out, we use the sum convention HijFj := ■ HijFj for 
spatial tensor indices. Restricting ourselves to dilute 
suspensions of slender dumbbells, we adopt the Oseen 
approximation for the hydrodynamic interaction tensor, 



i.e. 



ij 
H ij 



H. 



(as)(pr) 
ij 



(6nfJta) x 5ij, 
\X as -X pr \ 

(xr - 



N?N° 



(Al6a) 
(A16b) 



Su + 



xn(x 



as 
3 



x?) 



\x c 



(Al6c) 



where k = (8717/) 1 . To obtain from Eqs. (A15) closed 



stroke-averaged equations for {R 7 , N a }, we must 

a. perform a far-field expansion of the hydrodynamic 
interaction tensor; 

b. determine the internal forces F as , required to 
maintain the dumbbells' prescribed shape L a (t)\ 

c. expand the resulting equations in powers of 
and average over a stroke period [t, t + T]. 



a. Far-field expansion 



The Oseen tensor components given in Eq. (A16) 
are functions of the sphere separation vectors X as — X pr . 
By means of Eq. (A14c), we may decompose 

X as _ x pr = D ap + y (as){pr) ^ 



where similar to Eqs. (A12b) we have defined 

D ap := R a R p , (Al8a) 



Y (.° s )(pr) .- [(-l) s N a L a - {-l) r N p L p } 



(Al8b) 



Then, for a 7^ p, the Oseen tensor components (Al6c 
take the form 



Hi 



\D + Y\ 



Di + Yj Dj+Yj 
\D + Y\ \D + Y\ 



(A19) 



For clarity, we dropped superscripts here using the ab- 
breviations Y := Y {as){pr) and D := D ap . In the dilute 
limit, corresponding to \Y\ <C \D\ we may perform a far- 
field (Taylor) expansion of the tensor components ify. 
For this purpose we define 

H% := Hj (Y = 0) = ^ (Sa + AA) . (A20a) 
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where 



D, := 



A 

\D\ 



(A20b) 



is the unit vector in the direction of D. Reinstating up- 
per indices, the formal Tay 
Y = can be expressed as 



per indices, the formal Taylor expansion of JJ^J at 



Y 



/ j ij.k q ...k\ ki 
9=0 



(crs)(pr) 



■■Y, 



(as)(pr) 



where 



n i 3 ,k q 



.fci 



-y<9 fel ■■■d kq H° J 



(A21a) 



(A21b) 



D=D ap 



and dk '■= d/dD k . Explicit expressions for the expan- 
sion coefficients Hij,k q ...k^ with q = 1,2,3,4 are summa- 



rized in B The expansion ( A21 ) will be used in the next 



part to compute the interaction forces F as , and, later 
on, it will also be inserted into the exact equations of 



motion ( A15 1 



b. Internal forces in two-particle approximation 



We wish to determine the internal forces F as in 



Eq. ( A15 1 by means of an iterative procedure, restrict- 



ing ourselves to two-body interactions and assuming, as 
usual, that individual dumbbell swimmers are both force- 
free and torque- free, i.e., 



s 

= T[(y) := £ e ijk (Xp y 3 )F^ s , (A22b) 



where y = (t/ 7 ) is an arbitrary reference point. Substi- 
tuting Eq. dA22ab into Eq. \A22H we find that 



= e ijk (Xf - Xf)F% \ 



or equivalently 



= eijkNfFZ 1 . 



This implies that F as must be of the form 



r 



-r 



(A23) 



(A24) 



It thus remains to express the N unknown functions f al 
in terms of {R 7 , N a \. 

Shape- constraints.- To determine the unknown func- 
tions f 1 , we exploit the N independent shape con- 
straints 



L° 



(A25) 



Inserting the equations of motion for X as , we find the 
explicit condition 



i CT = E^° 



'(a2){pr) „(<7l)(pr) 



NPf pr . (A26) 



p , r 



Introducing the convenient abbreviation 

h {<js)(pr) ._ N a jj{crs){pr) jy-p ^ 

we can write Eq. ( |A26 1 as 



(A27) 



L° 



i \h {a2){,7r) - h {ai){,7r) 

r 

^2 \h { - a2){pr) - h {ai){pr) 



r 



p^cr.r 



epr 



(A28) 



Here, we have separated interactions within the dumb- 
bell a from those with other swimmers p ^ a. Using the 
force-free constraint (A24), Eq. ( |A28[ ) takes the form 



(A29a) 



with coefficient functions 



b ap := h {a2){pl) + h {al){p2) - [h {,j2) ( p2) + h {al){pl) ] 

(A29b) 



The N linear equations (A29a) determine the N un- 
known functions f pl by means of an iterative procedure. 



Iteration scheme.- Rewriting Eq. ( A29a) in the form 

(A30) 



T a yjp 

t*i = L V — fp l 



b aa t-^ b° 

we obtain the following recursive sequence 



- V — f pl 

Ucra J (n-1) 



^ b 

p^a 



(A31) 



Starting from the initial condition f p ^ = 0, the first it- 
eration gives the force generated by an isolated, shape- 
driven dumbbell 



fa. 

J 1 



(1) - b o 



(A32) 



The second iteration yields a correction due to pair in- 
teractions with other dumbbells, 



1(2) 



E 



f<x) ~ E 



p^a 



b°P ^ 

b aa bPP 

b ap LP 
b aa bPP ' 



(A33) 
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Similarly, one obtains from the third iteration 



We define 



'(3) 



U 7 
b aa 



pjtcr 



/(2)+EE 



LP 
bPP 

b a P bP y 



E 



bPP b vv 



( 2 ) 1 /-^i Z__/ 50-0- Jjpp ty>v ' 

p/o- V ^p 



(A34) 



The last term can be interpreted as a three-particle in- 
teraction correction. Let us assume that the system con- 
tains £7 = 1, . . . , iV dumbbells. Then, as evident from the 



'exclusive' summation in Eq. (A34), the iteration will 



approach a fixed point after N iterations, 



J(N+l) 



ral 



(A35) 



The fixed point frm corresponds to the exact solution, 
i.e., /(jv) is the internal force generated by a dumbbell 
in order to maintain its prescribed shape in the pres- 
ence hydrodynamic forces of JV — 1 other dumbbells. In 
the remainder, we shall restrict ourselves to considering 
one-particle and two-particle interactions, corresponding 
to m and f$ 



(1) ^ J(2)- 

Coefficients b ap ' . 



We still need to determine the co- 



efficients b ap from ( A29b ) . The 'diagonal' coefficients b° 



can be calculated exactly by noting that 



h {ol)(<72) 



4k 
3a 



= /j(^2)(al) = 



2k 



We thus have 



4k 



2L a 
3 a 



(A36a) 
(A36b) 

(A37) 



In order to determine the coefficients b ap with p ^ a, we 
need to use the far-field expansion (A21|. Defining the 
contraction 



allows us to write 



<J=0 



(A38) 



(A39a) 



where 



N <rap 

k\k 2 k 3 



-KTCrcrpp 
k 1 k 2 k 3 k 4 



nraacrp 

k 1 k 2 k 3 k 4 



N Zi N k 2 N k 3 + N k N k 3 N k 2 



(A40a) 



(A40b) 



KK 
KK 



NZ n N^N£NZ 3 +NZ 3 N^ 



K K , 

fel k 2 ' 



(A40c) 



KKKK + 
KKKK + 



KKKK + 



KK 



(A40d) 



With these abbreviations we find 



K p = 0, 

b7 = 0, 

b7 = h% ki L°L»N° k ? k2 , 

b%" = 0, 

4 a k^k^k^ki L /C1/C2&3&4. 



(A41a) 
(A41b) 
(A41c) 
(A41d) 



L p L p L p L a N ppp \ , 1 , 

k\k 2 k 3 ki\ ' 



(A41e) 



which can be used in ( A33 1 . 



c. Stroke- averaging 



Translational motion.- Inserting the ansatz (A24) 



into Eq. (A15 



is determined 



the motion of the position coordinate 



>y 



k i= lJ2 H t )(pr) f pr K- 



(A42) 



s,p,r 



It is convenient to consider 'internal' and external con- 
tributions separately by writing 



where 



p±o 



ja ._ 1 ^ ^(as)(ar) jar ^ 



J! 



(ers)(pr) j-prpjP 



(A43a) 

(A43b) 
(A43c) 



b ap 
1 



ft 



h ap 

kq . . . k\ 

(<xl)(p2) 

fel 



■■■Y, 



(<Tl)( P 2) 



(al)(pl) 



fel 



• Y, 



(<Tl)(pl) 



Y 



(<x2)(pl) 



fel 



(<r2)(p2) 



fei 



(ff2)(pl) 



(<x2)(p2) 



]}• 

(A39b) 



Here we have used that 



Using the force- free constraint and Eq. ( A36 1 , we find 
I a = 0, (A44) 
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i.e., the only contribution to the translation of swimmer a 
comes from interactions with the other dumbbells p 7^ a. 
Hence, we still need to determine the second contribution 
J° p from Eq. ( A43c I , which can be written in the form 



7°"P _ t!<rp njP fP 1 
°i "ij 3 ' 



(A45a) 



where 



1 



E ap := ~\H-- 
ij 2 3 

1 r /7 -(<r2)(pl) 

2 I ij 



(<rl)01) _ ff ("l)(p2) 



H 



n J 

(ff2)(p2)-| 



(A45b) 



Inserting the far-field expansion for the Oseen tensor, we 
obtain 



00 

"ij / j ij,k q ...k 1 ki...k q 
q=0 

up 



(Y) (A46a) 



with polynomials P k p k given by 



pop (y) •= ly( fTl )(^ )1 ) . 

ki...k q \ ' ' 2 


.y fc M)(pl) _ 




l v -(al)(p2) 

2 fc i 


y ^(<Tl)( P 2) + 




l v (tr2)(pl) 

2 fc i 


. y>2)(pl) _ 




l v (<72)(p2) 

2 fe i 


. r (^2)(p2) 


(A46b) 



where D := D ap := R a - R p , D := D ap /\D lTp \ and 



s = Dj P N? , r = D° P N P , q = N° N p 

denote the three possible pairwise projections of the 
relevant unit vectors N a , N p , and D P . Inserting 
Eqs. (A49l into (A48l yields the expression for J CTp that 



is given in Eq. (9c 



Change of orientation.- The exact equations of mo- 
tion for the orientation vectors N a read 



up 
k ' 



(A50a) 



where 



„(<x2)(p&) _ Tr{ul)(pb) 
n kj n kj 



(A50b) 



Using the force-free constraint (A22a), one obtains ex- 
plicitly 



G a k p = N p [-H^/ Kpl) +H ( k f ip2) + 

(u2)(pl) _ H (a2){p2) ] 



H. 



(A51) 



kj kj J £. 

Inserting the expansion for hydrodynamic tensor gives 



In particular, for q = 0we have P 



up 



p<?P 
r k lk2 

p^P 

k\k 2 k 3 



pop 
ki . . ./C4 



L p N p , 

fci ' 



and for q > 1 

(A47a) 
(A47b) 



[L*L*IS NlZks + 
L p L p L p N p N p N P 3 ], 



(A47c) 
(A47d) 



Since Eq. (A43a) already contains a sum over p, ne- 



glecting three-body effects means that, in order to com- 
pute Jf p , we should use f pl ~ f p { \ = LP/bPP in 



Eq. (A45a). After averaging (A45a) over period, we ob- 



tain at leading order of (i/\D 



H (ul)(pl) + Hyd iol)(p2) + H (u2){pl) _ H (p2)(p2) 



ETTUP 
ij,k q ...ki 

q=Q 

[-Y k (alKpV 



■■■Y, 



(ol)(pl) , Y (ul)(p2) 



■■■Y, 



(ul)(p2) 



Y, 



( CT 2)(pl) 



fei 



Y, 



( CT 2)(pl) _ Y (u2)(p2) 



fel 



■ • Y, 



(u2)(p2)l 



The polynomial terms in brackets are exactly those en- 
countered earlier in Eq. ( A39b |. Hence, the first two non- 



vanishing contributions come from q — 2 and q — 4. Ne- 
glecting three-body effects means that, similar to above, 
we should use f pl ~ f p ^ = L p /b pp . Hence, truncating 
after q = 4 we have 



1 I J <y\2 J p f p 

JOP ^ ± JSJP TTUP JSjUp \U ) ±J U 

i — 4 3 11 ij,k 3 k 2 k 1 ly k 1 k 2 k 3 ^ pp ■ 

where to leading order in A 



(L°) 2 LPLP 



bPP 8k 
The contraction is obtained as 



3^ 2 A 2 . . 
uia — sinf^^ 



(A48) 



(A49a) 



\N p HZNZ = ~^{N? {2s + Aqr - 10,r 2 ) + 

A(l + 2q 2 - 5s 2 - 5r 2 - 20qsr + 35s 2 r 2 )}. 

(A49b) 



L P L P 

- - n i3Mki I "k 1 k 2 b pp - n ij,fc 4 fe3fe2fcl 

1 u 2 a 3 L p 

4 [ N kik 2 k 3 ki ( L<J ) LP + N kik 2 k 3 k 4 ( LP ) ] ipp- 

Averaging this expression over a period, we find 



1 I T a\2 T p f p 

U i ~ ^3 It ijMk 3 k 2 k 1 ly k 1 k 2 k 3 ki b pp ■ 



The time average on the rhs. is the same as in (A49a| 



Exploiting the symmetry of lower indices of NZ^Tu, we 
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obtain 



4 ^v ■ nij,knlm^ m t n k ~ 4|£)|5 



A^(-3 + 6g 2 + 15s 2 + 15r 2 

-105sV + 60srq) + 
5A(3s + 6rq + 6sq 2 - 7s 3 - 21s?' 



-42<js 2 r + 63sV) 



Contracting with the orthogonal projector (6^ — N£ N[), 
see Eq. ( |A50a ), eliminates terms proportional to N? , 
thus yielding Eq. (9b I. 



Appendix B: Partial derivatives of the Oseen tensor 



Second order derivatives.- The second order deriva- 
tives, normalized by n\ with n — 2, are defined by 



and read explicitly 



H, 



ij,kn 



2!|D| 3 

+3 x ( 



D n DkSij — D n Dj8ik — D n DiSjk — 

DiDj5 nk - DkDjSni - DiD k 5 nj ) 

+3 x 5 x Dnbkb.Dj] . (B5) 



This part summarizes the partial derivatives of the 
Oseen tensor that are required in the derivation of the 
far-field, stroke-averaged equations of motion Q, see 



Eq. ( A21 1 in A2 



Consider the distance vector D = (Dk), its associated 
unit vector (D k ) and orthogonal projector (11^), given 
by 



-Dfc 
\D\ 



n, 



& •= 



DjDi. 



(Bl) 



We wish to compute the partial derivatives of the Oseen 
tensor 



Hi 4 '■ = 



D\ V 



1 + D.Dj 



where k := (87r/i) 1 . Abbreviating d k 
have 



d k \D\ = = b k 



Sjk _ DkDi _ ITfe 
\D\ |D|3 ~ \D\ 



(B2) 
d/dDk, we 

(B3a) 
(B3b) 
(B3c) 



First order derivatives.- A straightforward calcula- 
tion gives 



H, 



ij,k 



dkH lJ 
D 



(-DkSij + Dj5 lk + DiSjk 



|£>l 2 



(B4) 



Third order derivatives.- Similarly we find for the 
third order derivatives 



Hij,kni '■— 7ydid n dkHij 



the explicit representation 



H, 



ij,knl 



3!|D| 4 

3 x [Di(S n kSij - 5 nj Sik - 5 ni 5 jk ) + 
D n (SikSij - SijSik - SuSjk) + 



Dk(Sl n Si 



Sn6. 



Ij Oni 



SliSnj) 



Di(5l n Sjk + 0~lj5nk + SlkS n j) - 
bj(6ikS ni + Si n Sik + 0~liO~ n k)] 

+3 x 5 x ( 

-btbnbkSij + bibnbjdik - 
bib„bi8 jk + bib.bjSnk + 
bib k bjS m + bibib k s nj + 
b k bib 3 s ln + bn^bjSik + 
b n b k b j s H + b n b k b i 5 lj ) 



-3 x 5 x 7DiD n D k DiD,j 



(B6) 



Fourth order partial derivatives.- Finally, the fourth 
order derivatives, defined by 



H, 



ij,knlm 



dmdldndkHij 
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are obtained as 



K f r 

= 4wr x[ 



Sml{S nk Sij SnjSik SniSjk) 
+8mn(filkfii 3 - SljSik - SuSjk) 

+5 m k{3in$i 3 - 5i 3 5 ni - 5 U S nJ ) 

— 5 m i(5l n 5jk + SijS n k + SlkS n j) 

— 6 m j(5lkS n i + Si n 5ik + SuSnk)] 

+ 3 x 5 x [ 

D m Di(5i n 5jk + Sij5 n k + SlkS n j) 

+D m Dj(Si k 5 ni + Si n S ik + SiiS n k) 

+DiDi(S mn 8jk + SmjSnk + S m kS n j) 
+DiD J (S mn S i k + S mt S n k + S m kS n i) 
+D n Di(5 m lSjk + SmjSlk + 8mk$lj) 
+D n Dj(S m lSik + Smifilk + SmkSli) 

+DkDi(5 m i5j n + S m jSi n + S mn Sij) 
+D k Dj(S m i6 

+DiDj(S m iSkn + 5 m k5in + S mn 8ik) 
—D m Di(5 n k8ij — 5 n j5ik — S n iSjk) 
—D m D n (Sik5ij — SijSik — SuSjk) 
—D m D k (Si n 5ij — SijS ni — Sn5 n j) 

— D n Dk(S m l5ij — S m iSlj — 5liS m j) 
DlDk{SmnSij SmiSnj SniSmj) 

-DiD n (S mk 5 l: j - S mi S jk - S ik 5 mj )] 
+ 3 x 5 x 7 x [ 

D m DiD n D k 5 l0 — D m DiD n DjS ik 
—DmDiDnDiSjk — D m DiDiD 3 S nk 
—D m DiD k DjS ni — D m DiDiD k S n j 
—D m D k DiDj5i n — D m D n DiDjSi k 
-D m D n D k DjSu - D m D n D k Di5ij 
—8 m iD n D k DiDj — S mn DiD k DiDj 
—5 mk DiD n DiDj — S m iDiD n D k Dj 
-S mj DiD n D k Di] 
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